More investigation into alpha cost
DRAFT / WORK IN PROGRESS
Josh Nugent
# Load required packages and data
library(tidyverse)
library(pander)
library(fitdistrplus)
library(AER)
dat <- readRDS("t_test_data_sbs_pt01_10sub_10clust_10000sims.rds")
pval_b1_model <- unlist(dat[1,])
est_sb2 <- unlist(dat[2,]) # let's assume anything less than 1e-6 is the weird ones
tiny_est <- est_sb2 < 1e-6
resid_var_model <- unlist(dat[3,])
MSB <- unlist(dat[5,])
MSE <- unlist(dat[6,])
clustvar0 <- unlist(dat[7,])
clustvar1 <- unlist(dat[8,])
ftest_pval <- unlist(dat[9,])
pval_tt_eq_var <- unlist(dat[13,])
pval_tt_uneq_var <- unlist(dat[14,])
df_tt_uneq_var <- unlist(dat[15,])
beta0 <- unlist(dat[17,])
beta0_se <- unlist(dat[18,])
beta1 <- unlist(dat[19,])
beta1_se <- unlist(dat[20,])
plotdat <- cbind.data.frame(tiny_est, pval_b1_model, est_sb2, resid_var_model, #modelNLME,
MSB, MSE,
clustvar0, clustvar1, ftest_pval,# ttest_eq_var, ttest_uneq_var,
#clustmeans,
pval_tt_eq_var, pval_tt_uneq_var, df_tt_uneq_var,
beta0, beta0_se, beta1, beta1_se)#, test)
plotdat$icc <- MSB / (MSB + MSE)Data-generating mechanism
- Two arms
- 10 clusters per arm
- 10 subjects per cluster
- \(\sigma_b^2 = .01\) in both arms
- \(\sigma^2 = 1\) in both arms
- \(\beta = 0\) in both arms (no treatment effect)
Gathering data on how a ‘naive’ NLME model estimates those parameters, as well as what equal- and unequal-variance t-tests show.
When we get the tiny \(\sigma_b^2\) estimates, ratio of MSB to MSE from original dataset is small
Mean square within clusters (clusters \(i = 1, 2, ..., k\); subjects \(j = 1, 2, ..., n\) in each cluster): \[ MSE = \frac{1}{k(n-1)} \sum_i \sum_j (y_{ij} - \bar{y}_{i.} )^2 \] Mean square between clusters: \[ MSB = \frac{1}{k-1} \sum_i \sum_j (\bar{y}_{i.} - \bar{y}_{..} )^2 \]
I looked at a lot of relationships between variables, and the ratio \(MSB / MSE\) was the one that popped out as a clear predictor of a tiny estimate.
cutoff_MSBMSE <- min(MSB[tiny_est==F] / MSE[tiny_est==F])
cutoff_ICC <- min(plotdat$icc[tiny_est==F])
ggplot(data = plotdat) +
geom_histogram(aes(x = MSB/MSE), color = "black", fill = "grey") +
labs(title = "Ratio around .1 highly determinative of a tiny estimate") +
facet_grid(rows = vars(tiny_est), labeller = label_both)
ggplot(data = plotdat) +
geom_histogram(aes(x = MSB), color = "black", fill = "grey") +
labs(title = "MSB alone is less predictive") +
facet_grid(rows = vars(tiny_est), labeller = label_both)
ggplot(data = plotdat) +
geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) +
labs(title = "Ratio around .1 highly determinative of a tiny estimate")
ggplot(data = plotdat) +
geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_abline(aes(slope = 1/cutoff_MSBMSE, intercept = 0)) +
labs(title = "More variability in that finding among the tiny estimates") +
facet_grid(rows = vars(tiny_est), labeller = label_both)
ggplot(data = plotdat) +
labs(title = "No change in the pattern as MSE varies") +
geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))
ggplot(data = plotdat) +
labs(title = "No change in the pattern as MSB varies") +
geom_point(aes(x = MSB, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))
ggplot(data = plotdat) +
geom_point(aes(x = clustvar0/clustvar1, y = MSB/MSE, color = tiny_est), alpha = .3, stroke = 0) + xlim(c(0,4))
labs(title = "Pattern seems invariant to ratio of between-cluster variances")## $title
## [1] "Pattern seems invariant to ratio of between-cluster variances"
##
## attr(,"class")
## [1] "labels"ggplot(data = plotdat) +
geom_point(aes(y = MSB/MSE, x = clustvar0 - clustvar1, color = tiny_est), alpha = .3, stroke = 0) +
labs(title = "Pattern seems invariant to difference in between-cluster variances")
ggplot(data = plotdat) +
geom_point(aes(x = MSE, y = icc, color = tiny_est), alpha = .2, stroke = 0) +
labs(title = "Using the ICC instead of the ratio might also work, results are the same") +
geom_hline(aes(yintercept = cutoff_ICC)) +
facet_grid(rows = vars(tiny_est), labeller = label_both)
ggplot(data = plotdat) +
geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) +
labs(title = "Using the ratio") +
geom_hline(aes(yintercept = cutoff_MSBMSE)) +
facet_grid(rows = vars(tiny_est), labeller = label_both)
sum((MSB[tiny_est==T] / MSE[tiny_est==T]) > cutoff_MSBMSE)## [1] 628sum((plotdat$icc[tiny_est==T]) > cutoff_ICC)## [1] 628ggplot(data = plotdat) +
geom_histogram(aes(x = pval_tt_uneq_var), color = "black", fill = "grey") +
facet_grid(rows = vars(tiny_est), labeller = label_both) +
xlab("p-value for welch t-test, allows unequal variances (results the same with eq vars)") +
labs(title = "Small p-values on a t-test correlated with tiny estimates; makes sense, more power")
What to make of this? Just a signal-to-noise issue?
What about the cases that violate that boundary line?
First, find some wild outliers and examine them…
library(ggplot2)
library(pander)
outliers <- plotdat[ MSB/MSE > cutoff_MSBMSE+.05 & tiny_est == T, ]
outindx <- c(779, 4719, 6551, 7571)
plotdat$outlier <- FALSE
plotdat$outlier[outindx] <- TRUE
ggplot(data = plotdat) +
geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est, stroke = outlier, alpha = outlier)) +
labs(title = "Looking at these four most extreme outliers...") +
geom_hline(aes(yintercept = cutoff_MSBMSE))## Warning: Using alpha for a discrete variable is not advised.
pander(outliers[,c(4:7,13,15)]) # Clue!| resid_var_model | MSB | MSE | clustvar0 | beta0 | beta1 | |
|---|---|---|---|---|---|---|
| 779 | 0.9266 | 0.1417 | 0.9401 | 0.04995 | -0.3351 | 0.5034 |
| 4719 | 0.9447 | 0.1444 | 0.9594 | 0.03254 | 0.1923 | -0.5114 |
| 6551 | 0.9413 | 0.1563 | 0.9545 | 0.06521 | -0.3572 | 0.5499 |
| 7571 | 0.8346 | 0.1555 | 0.856 | 0.05251 | 0.4466 | -0.6063 |
ggplot(data = plotdat[tiny_est==T,]) +
geom_point(aes(x = MSE, y = MSB/MSE, color = abs(beta0 - beta1)), stroke = 0, alpha = .5) +
labs(title = "Outliers more common with large differences in model beta coefs...") +
geom_hline(aes(yintercept = cutoff_MSBMSE))
ggplot(data = plotdat[tiny_est==T,]) +
geom_point(aes(x = MSE, y = MSB/MSE, stroke = outlier, color = abs(beta0 - beta1) < .45), alpha = .5) +
labs(title = "... easier to see when discretized.") +
geom_hline(aes(yintercept = cutoff_MSBMSE))
Still, outliers not fully explained.
Other looks at the outlier data sets haven’t turned up any more clues (yet).
# Something with the underlying data sets?
outlier_datasets <- dat[16,outindx]
ggplot(data = outlier_datasets[[1]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[2]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[3]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))
ggplot(data = outlier_datasets[[4]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))